function [sample_means sample_covars] = computeSufficientStatisticsUncertainty(data, vars, covar_type, reg)


K = size(vars.Z, 2);
[N D] = size(data.Y);

sample_means = zeros(K, D);
sample_covars = zeros(D, D, K);
for k = 1:K
    w = vars.Z(:, k);
    mydata.Y = vars.R(:, :, k);
    mydata.C = vars.E(:, :, :, k);
    mydata.YYt = vars.RRt(:, :, :, k);

    [mu sigma] = computeMLEGaussianUncertainty(mydata, w, covar_type, reg);
    sample_means(k, :) = mu;
    sample_covars(:, :, k) = sigma;
end
